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^ ■ We study the electrical transport of a harmonically-bound, single-molecule Ceoshuttle operating 

^ ' in the Coulomb blockade regime, i.e. single electron shuttling. In particular we examine the 

, dependance of the tunnel current on an ultra-small stationary force exerted on the shuttle. As 

an example we consider the force exerted on an endohedral NQCeoby the magnetic field gradient 
generated by a nearby nanomagnet. We derive a Hamiltonian for the full shuttle system which 
includes the metallic contacts, the spatially dependent tunnel couplings to the shuttle, the electronic 
and motional degrees of freedom of the shuttle itself and a coupling of the shuttle's motion to 
a phonon bath. We analyse the resulting quantum master equation and find that, due to the 
exponential dependance of the tunnel probability on the shuttle-contact separation, the current is 
highly sensitive to very small forces. In particular we predict that the spin state of the endohedral 
S-H ' electrons of NQCeofn a large magnetic gradient field can be distinguished from the resulting current 

signals within a few tens of nanoseconds. This effect could prove useful for the detection of the 
endohedral spin-state of individual paramagnetic molecules such as N@C6oand P@C6o, or the 
detection of very small static forces acting on a Ceoshuttle. 
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PACS numbers: 72.70.+m,73.23.-b,73.63.Kv,62.25.+g,61.46.+w,42.50.Lc 



I. INTRODUCTION 



The determination of the spin-state of a single electron spin has attracted much attention recently Q . Following the 
ground breaking experiments of Park et al., in the observation of quantized motional transitions in the conductance 
of a single-molecule Ceo transistor interest in quantum electromechanical systems, or QEMS, has blossomed 0- In 
this work we wish to examine whether a QEMS, with the inclusion of spin degrees of freedom, or a Spin-QEMS, can 
^nI" ■ be used to detect the spin-state of the included spin. We have in mind the particular case of Nitrogen or Phosphorous 
endohedrally doped Ceo, N@C6oor P@Ceo- This material has previously been shown to possess unpaired endohedrally 
trapped electrons, in a quartet ground state, S = 3/2, and which display exceedingly long transverse relaxation (or 
' T2), times. This material has been considered for use in quantum information processing devices but may have uses 
in other areas related to classical spintronic memories. 

As a first exploration one can consider the original experiment of Park et al. using the paramagnetic endohedral, like 
PQCeo, instead of Cgo- It is known that the electronic spin degrees of freedom reside in the endohedral electrons which 
are concentrated isotropically at the centre of the Cgocagc in P@C6o ;4|. In the experiment of 0], the Ceooscillates in 
' O I a Harmonic potential at THz frequencies. For the Spin-QEMS system of a P@C6omolecule as a movable island in a 
single molecule transistor, one must be capable of engineering a coupling between the state of the endohedral spin and 
the conductance properties of the transistor. This can be accomplished through the imposition of a large magnetic 
field gradient dB/dx, along the plane of the transistor. The interaction between the endohedral spin and this gradient 
field will impart a small force on the molecule which will in turn, cause a slight shift in the equilibrium position within 
the damped Harmonic oscillator trapping potential, shifting the molecule towards one of the contacts. This small 
jj] ■ shift will lead to a change in conductance properties. It will be shown that if the configuration is asymmetric, i.e. the 
mobile island is closer to one lead than the other, the resulting current depends exponentially, on this small spatial 
displacement. In the case of a tight trapping, observed in Park et al., Q, the resulting spin-dependent change in 
conductance is extremely small and unobservable. If we instead configure for the largest achievable magnetic field 
gradient (obtained for instance in MRFM devices), by placing a nanoscopic permanent magnet near to the mobile 
island, we also find that due to the large difference between the energy scales of the endohedral electronic spin flip 
(GHz), and the island's harmonic motion in 0, (THz), the island's spin-dependent position shift is minute and is 
unobservable in the current. 
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Guided by these observations we instead consider the more "floppy" Spin-QEMS of a shuttle device where the 
restoring oscillator now operates in the GHz frequency regime (see Figure MRFM cantilevers fabricated out of 
single Si crystals have recently been generated which have resonant frequencies in the GHz range ||- The resulting 
device could be similar to that outlined by Isacsson [fj. It may also be possible to engineer non-conducting flexible 
linker molecules which will connect the shuttle molecule to the metallic contacts and have resonant frequencies in the 
GHz range. 

We thus consider the system of a movable shuttle or island, which is harmonically bound to oscillate between two 
metallic contacts in the presence of significant motional damping. We further assume that we operate in the Coulomb 
blockade region where only one electron can reside on the shuttle at one time. For the particular case of an endohedral 
N{P}@Cgoisland, as in y|, we make the assumption that this injected electron distributes itself isotropically around 
the surface of the Cgocage thus averaging out any magnetic dipole coupling between the spin state of the injected 
electron and the spin state of the endohedral electrons. We also assume that the injected electrons are not spin 
polarised and the resulting force experienced by the shuttle due to the interaction between the spin of the injected 
electrons and the magnetic gradient will average out while the force on the shuttle due to the endohedral electrons 
interacting with the magnetic field remains static. We therefore neglect the random spin-force contribution of the 
injected electrons in what follows. 

Before going into the model in detail, we can understand, in a simple but surprisingly robust manner the underlying 
reasons why we expect to achieve an exponential dependence of the current on the island's displacement from its 
equilibrium position. The current through a simple source-fixed island-drain, two-stage sequential tunneling process 
in terms of the left and right tunnel rates 7£, jr, is given by 

1L1R ,^ 



1R + 7L 



If we now take into account the exponential dependence of these tunnel rates on the position of the (now mobile) 
island, and letting S (positive) be the small displacements of the island from its equilibrium position to the right (to 
the drain), we can set 

lL-lle~ U , l R ~l R e+™ , (2) 

and thus the right tunnel rate, from the island to the drain, increases due to the smaller separation while the left 
tunnel rate decreases due to the larger separation. This is reversed if the motion is to the left (S negative), or towards 
the source: 

7i ~ 7 °e+W 7i? ~ 7 °e-^l . (3) 

Inserting these into Q, we easily get: 

■positive displacement -.0~ — A|<5| i ^,0 r -t-A|6 i i l? p +2\\8\ 

^ _ ° _ /L c ~T IR 1 - _ 1 ~r r L 

^negative displacement ,y0g+A|(5| _|_ ^y0 g — A|<5 g+2A|<5 _|_ JP 1 ^ ' 

where F = 7^/7° . If F = 1 then fi = 1, and there is no difference between the currents in the two cases. When 
F e +2A l <5 l, ~ e +2A|<5|^ w \j[\ e if p ^ e~ 2A l' 5 l, f2 ~ e _2A ' (5 ', and in cither of these asymmetric cases the two currents 
can be very different from each other if \\5\ is appreciable. 

In the sections below, we analyse the classical and quantum models of this system. Following the above sim- 
ple rationale, wc can find a realistic parameter range where the right-hand tunnel current shows a large (exponential), 
dependance on small spin-dependent island position displacements. This is confirmed below using low-resolution 
numerical simulations of both the "semiclassical" and quantum systems, and through high-resolution quantum 
simulations. Besides demonstrating the two different current "signals" corresponding to the spin force, we must also 
determine the noise present within the system to estimate the signal to noise ratio or so-called Fano factor. Following 
standard methods wc derive the current spectral noise density and from this the signal to noise properties of the 
coupled system. Using all this we finally deduce the measurement time required to distinguish the spin state from 
the current signals and can find device parameters where this acquisition time will be several tens of nanoseconds. 



II. MAGNETIC FIELD GRADIENT 



We now estimate the magnitude of the maximum magnetic field gradient wc can achieve at the location of the 
shuttle. We follow a similar derivation in Magnetic Resonance Force Microscopy (MRFM) 0. The Hamiltonian for 
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the shuttle's cndohedral spin, S, interacting with a magnetic field which varies along the x— axis (see Figurc[3J|, B(x), 
is H = —fix ■ B = -~HB9eB{x)S z , while the resulting force on the shuttle is given by F spin = pLgg e S z dB(x) / dx. If we 
consider placing a nanoscopic ferromagnetic spherical particle of radius R a distance d from the shuttle, one has 

- _ 3n(m ■ ft) - rh 

JDFerro — , (O) 

An r A 

where fio = Air x 10 _7 H/m, rfi is the magnetic moment of the ferromagnetic particle, r is the vector connecting the 
particle to the shuttle and r = \ f\ = R + d. For a sphere the magnitude of the magnetic moment is \m\ = A/3irR 3 M , 
where M is the magnetization of the ferromagnetic material. Following this we have Bfctto = — fioR 3 M/r 3 z. As the 
shuttle moves along the x— axis, the separations d and r, between the shuttle and the magnet change. Setting the 
distance between the shuttle and the magnet's centre to be r = R + do + x = x + x, where x is the small oscillation 
of the shuttle, do is the equilibrium separation, and x = R + do, one can derive 

^ m ms z {rV 

x \ x J 

For an iron ferromagnet, [IqM = 2.2Tesla, and choosing c? = 5nm, R = 5nm, x = lOnm, S z = ±3/2 (for P@Cgo), 
one obtains F spin ~ ±1 x 10~ 15 N. 

Using this we can re-examine the spin Hamiltonian and expand to first order in x about the equilibrium shuttle 
position, H = —fiBg e S z B z (r) = —fiBg e S z B z (x + x) ~ Ho — F sp i n x + 0{x 2 ), where Ho is a constant which we drop. 
From above we set 

B spin \^z^ : (7) 

where 



3^j, e g e HoM ( R 



(8) 

x \ x J 

The Hamiltonian for the Coloumb-energy of the charged island/shuttle can be written as 

H shuttle electronic — flldlC C , (9) 

where c represents the operator which annihilates an electron on the shuttle, and c^c is the number of injected electrons 
on the shuttle, which in the blockade situation is cither or 1. 



III. ELECTROSTATIC FORCE 

In most studies of electron shuttles 0, one considers the shuttle to be driven primarily by the electrostatic force 
exerted on the shuttle through the interaction of the shuttle charge and an electric field present between the source 
and drain metallic contacts created by a potential difference V = Vl — Vr, across the inter-contact gap D. This gives 
an electric field of strength E = V/D. 

In a standard treatment of electron shuttling the shuttle acquires the electron at the average displacement of 
zero slightly towards the source and continues to move towards the source. It then undergoes electrostatic force to 
be propelled back towards the drain and loses the electron at a slightly displaced position towards the drain while 
continuing to move towards it. Again under electrostatic force the shuttle moves back towards the source and repeats 
the process. In this case no harmonic restoring force is required to sustain the shuttle instability. However in our 
case, in the situation where our shuttle is operating in the Coulomb blockade regime, the shuttle can acquire a single 
excess negative charge near the source and thus be forced electrostatically to the drain, but once the shuttle dumps 
its electron at a displaced position near the drain there is no electrostatic force which can bring this shuttle back to 
the source contact to repeat the shuttle sequence. Thus in our case a Harmonic restoring force is required to sustain 
the shuttling action. 

In more detail, Fes = qE = -eE, where E = -dV/dx - (V R - V L )/D = -V/D, if V L = -V R = V/2. Setting 
the potential due to this electrostatic force to be <j)Es{x), Fes = —d(j)ES / dx, and 4>es{ x ) = — J Fssdx = —qEx = 
\e\Ex = -\e\Vx/D. Setting fj = V\e\/D, we have 



Hes = -fj x 



(10) 
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IV. HARMONIC OSCILLATOR 

The Hamiltonian of the Harmonic oscillator is dealt with in the standard manner. Setting the annihilation operator 
for motional quanta to be, s/2a = x/xq + ip/po, where x = wh/JmuJo) and p = V muj o^ are the extents in x and p, 
of the ground state wave function of the harmonic oscillator, we have 

„ P 2 . k 2 p 2 muj 2 2 t , . 

Hho = t, 1 tt x =7, \- -TT- x = niuo a> a . (11) 

2m 2 2m 2 

V. COMPLETE HAMILTONIAN 

We are now in a position to write out the complete Hamiltonian for the drain, source and motional, spin and 
electronic degrees of freedom of the shuttle: 

H = huJrJc (12) 

- X S z (^ + a)/V2 (13) 

- r){a) + a)Jc/V2 (14) 
+ Filuqo) a (15) 

+ ^E<4,X as) 

ki 
k r 

+ hJ2 "k b b\b kb + (]T a*b ki + h.c.) (18) 

kb ki 

+ ^(T4^ L (fi)o^ct+h.c) (19) 
+ ^(tE fi (i)«tct + h,.) (20) 

fctr 

where \ = X^o- an d T) = fjxo. In the above: (|12|l is the Coloumb self-energy of the shuttle, (|13fl is the spin-dependent 
potential energy, (|14|1 is the electrostatic potential energy, (jl5(l is the Harmonic oscillator self-energy, l|16(l and (|1 Tf> arc 
the self-energies of the Fermi baths in the source and drain contacts, (|18|) is the self-energy and coupling of a phonon 
bath to the motion of the shuttle, H19|l is the tunneling interaction between the source contact and the shuttle, while 
(|20|l is the tunneling interaction between the shuttle and the drain contact. 

We note the operators El/r(x), in the tunneling terms. The amplitude for tunneling depends exponentially on the 
spatial separation between the source/drain and the shuttle. One can assume the forms 

E L {£) = l)! 2 ^ XLX , E R (x) = ^ 2 e + ^ x , (21) 

where x is the displacement of the shuttle from its equilibrium position. In the following we assume the source and 
drain are identical in shape and material and thus Xl = A/j = A, where for Gold contacts, A -1 ~ 3 Angstroms. As 
above we set A = Xxo to get 

El/r = lif R cxp(TA(a t + a)/y/2) . (22) 

VI. DISPLACEMENT PICTURE 

As we mentioned above, the small static spin force exerted on the shuttle, in the absence of the electrostatic force, 
will cause a slight shift in the shuttle's equilibrium position within the Harmonic trapping potential. We now move 
to a frame of reference where we shift this displacement back to the origin. This can be done by considering the 
transformation 

p^p 1 = D(a)prf{a) , (23) 
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where D(a) = exp(aa^ — a* a), is the standard displacement operator for operators satisfying [a^ ,a] = 1. By choosing 
a = —(xS z )/(^/2huJo), which corresponds to the physical displacement of the shuttle's equilibrium position due to 
the spin force F spin , of 

6=-^- , (24) 
we can apply the displacement transformation (|23[) . to the above full Hamiltonian to get 



H' = hu I c^c±Sr 1 c i c (25) 

- J7(a f +o)c + c/V2 (26) 

+ Hu>o<x 'a (27) 

+ hJ2uh a h a k, (28) 

+ hJ2»l4>Z ( 2 9) 

+ ft£ Wfci 6^ + atfe fcl + h.c.) (30) 

A;;, ki 



+ Y,(T k L t E L (x')a L kti c^+h.c.) (31) 
+ £(T££ R (*0< r c t +h.c.) (32) 

fetr 

where now 

= 7l /2 exp (-A[(at + a)/V2 T <5]) = ll /2 E L , (33) 
^fl(^) = iT exp (+A[(a t + a)/V2 T S]) = 1r /2 E r , (34) 

where 6 = 5/xq. 



VII. MASTER EQUATION 

Assuming that the time scales for the couplings to the source, drain and phonon baths are faster than those 
associated with the Harmonic oscillator and electrostatic potential, one can derive the following master equation |l0j : 

^ = -iu [a1a,p]+il[(at+a)c<c/V2,p] (35) 

+ 1l {h(ui T 5v)V[c*E L ]p + (1 - T 5r t ))V[cE L ]p} (36) 

+ 1R {/b(wj T ^[c+EjJp + (1 - f R (uj T S V ))V[cE R ]p} (37) 
+ k(N + l)V[a]p + KNV[a ] ]p , (38) 

where the super-operator V is defined to be V[A]B = ABA* - \{A^ AB + B A) . The mean excitation of the thermal 
bath is N = 1/ exp(hujo/kBT — 1) and f(x) = (exp(x / (ksT)) + 1) , is the Fermi factor where 

f L (uiT&ri) = f(urT8Ti-{iL)=f(uiTtTi-V/2) , (39) 
fnfai T Srj) = f{m T Sr, - m ) = f( Ul T St] + V/2) (40) 

with a dependence on the bias voltage (through the chemical potential). Notice the dependency of the Fermi factor 
on the spin force. If the displacement factor 5 is large enough, the Fermi energies corresponding to the two spin 
directions will become significantly different and large enough to be observable at low temperatures. If one supposes 
that the bias voltage can be tuned to be between these two energies such that only one Fermi factor (corresponding 
to spin up say), is near unity while the other is nearly vanishing, it should be possible to differentiate the two spins 
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based on the detection or absence of current. However with such low bias, we believe the noise in the system will be 
overwhelming and we thus proceed by assuming that the system is in large bias such that this effect is negligible. 

Now at milli-Kelvin temperatures, fc^T ~ 0.2/icV, while u>i ~ 5meV. So if V > lOmeV then fi — > 1, while f R — > 0, 
giving the simpler form for the Master equation: 

^ = -tw [o t o J p]+<2[( o t+o)c t c/v^,p] (41) 

+ lL V[c ] E L \p + lR V[cE R ]p 
+ k(N + l)V[a]p + nNV[a f ]p . 

VIII. EXPECTATION VALUES 

From the master equation 1)410. we can derive the following differential equations: 



(42) 



^=-o(f)-f(f) , (43) 
at po 2 xq 

^ = -M-) + F E3 ^c)-^) , (44) 

dt Xq 2 po 



where 



j? - V \ e \ \ Fe ^\ (a* 
Fes = jr— = (45) 

Dpo Po 

These equations can be further tidied by going to natural time units associated with the Harmonic oscillator, i.e. 
t = uiot, and by defining 

Il/r k _ F ES (ar , 

1L/R = , « = 7, > XES = , (46) 

u> 2ui ^o 

we get 

d{ch) 



dr 

d(x/x ) p _ x 



re 



dr p x 

d(p/po x 



(48) 



= -(-)+Xes(Jc)-k{^-) . (49) 
dr xo po 

The above equations cannot be solved as they are not closed. One can break the correlations in the semiclassical 
regime of high damping. When this is done we finally obtain: 

*M^L = 7 ie ± 2 ^e~ 2A <^>(l - (etc)) - j R e^ 2XS e +2X ^/^(^c) , (50) 
dr 

d(x/x ) p _ x 

dr po x 

d(p/p ) . x _ f p 

—j = -< — > +Xes(c 1 c) - k(— ) . 

dr x po 

while the complete master equation in these rescaled parameters takes the form, 

^ = -i\<Ja, p] + ixE.s[(a f + a)Jc/V2, p]/h + V[c) ' E L ]p + V[cE R ]p + k(N + l)V[a]p + nNV[a f ]p . (51) 

UT 
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where 

S il/i = A /^I(c/ct)e X p[TA5±A(a + a t )/\/2] . (52) 

In what follows we take the phonon bath to be at zero temperature N = 0, and that the typical phonon bath excitation 
will have energies much smaller than that of the oscillator, i.e. below 1GHz . 

We will, in particular, be interested in the "semiclassical" and fully quantum steady-state behaviour of the right 
tunnel current, or the expected value of the operator: 

i R ee ^-i R ee 7 fle =F2A« e +2Ax/«o c t c = ^ R e^ e +2Xs Jc , (53) 
e 

where x = x/xq, p = p/po, and i Rl is the re-scaled average current flowing from the moving island to the drain. We 
wish to determine whether (i R (t — > oo)) = (i R ), is significantly different in the two cases when the endohcdral spin 
is up or down. However it is not enough that the two currents (endohcdral spin up/down), arc different but we must 
also then determine the steady-state (or DC or zero- frequency) , quantum noise associated with these two steady-state 
currents. This DC noise will allow us to determine how distinguishable the two current signals are. This then will 
determine the minimum length of signal acquisition time one can average over to achieve a signal to noise ratio which 
will allow the endohcdral spin states to be distinguished to a high level of confidence. In the following subsections 
(|IX F3I we discuss the initial "semi-classical" and quantum mechanical steady-state for the coupled electronic-vibronic 
shuttle system, paying particular regard to the change in right-tunnel-current with endohedral spin-state. This is 
done in a parameter range of the system where we are certain that our numerics are very precise and faithfully 
simulate the system. We find a parameter regime where (i R ) / ' (i R ) can be substantial, indicating that the endohedral 
spin-state measurement could be possible. In subsection IjlX C|) . wc expand the size of the numerical simulation to 
cover a much larger range of parameters (see parameter set (B) in Table I), which models a much larger electrostatic 
drive force xes- This is done using the method of quantum trajectories |llL Il2l ITflj. This type of simulation is very 
costly and only a small number of such simulations were performed. We found that when the electrostatic driving 
force is increased, the ratio (i R )/ (i R ) remains significant. Further the highly accurate quantum trajectory simulations 
corresponds very well with the lower-resolution steady-state results indicating that our steady-state numerics could 
be trusted to quite large values of the electrostatic force. In subsection (|IX D|) . we derive the analytical expressions 
for the noise spectra and numerically derive the DC-noise component over a limited (but precisely modeled), range of 
parameters. We reproduced the expected values for the DC noise of a standard, spatially fixed quantum dot, in the 
case when A = 0, i.e. no coupling between the electronic and vibration degrees of freedom. We further found a slight 
decrease in this noise with the amount of electrostatic driving and degree of the exponential coupling. Finally, in 
subsection (|1X K[) . we gather together all this information to estimate the measurement time required to distinguish 
between the current signals due to endohedral spin up (down). 



IX. NUMERICAL ANALYSIS 



A. Physical Parameters 



We now determine the typical values of the dimensionless parameters appearing in the system's dynamics. We will 
take, mc 60 = 1.2 x 10~ 24 kg, ujq — 2ir x 1GHz, the potential difference between the contacts to be V = lOmcV, and 
the spatial separation between the contacts D = lOnm. From this wc have x ~ 1.2 A, 5 ~ 0.18, giving the physical 
displacement S ~ .2 A. The dimensionless displacement of the harmonically bound shuttle due to the electrostatic force 
alone in units of xq is Ses ~ I-^esI /(wwqXo) ~ 28, while A ~ 0.4. From symmetry conditions one might suspect that 
we cannot gain any advantage from the exponential dependence of the tunneling amplitudes on the shuttle-contact 
separation if the shuttle's equilibrium position is exactly midway between the contacts (effectively when 7l = 7r). 
Thus we choose jl ~ 0.1, while j R ~ 0.9. We choose R ~ 2, thus setting the system close to critical damping. We 
also set N = 0, or assume a zero temperature vibronic bath, throughout. The example re-scaled parameter values 
that we study below are summarized in Table I. 
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A 


6 


Xes 


K 


1L 


1R 


A 


0.4 


0.18 


28 


2 


0.1 


0.9 


B 


0.4 


0.4 


2 


2 


0.1 


0.9 


C 


0.4 


0.4 


10 


2 


0.1 


0.8 



TABLE I: Parameter values for the model (A) typical physical parameters, (B) initial model parameters studied in "semi- 
classical" and quantum steady state, numerical analyses, (C) parameters used in "semiclassical" , low-resolution quantum 
steady-state, and high-resolution quantum dynamical (quantum trajectory), numerical analyses. 



B. Semi-Classical and Quantum Steady-State 

From Eq. (|50|l . one can obtain the steady-state solution for the expectation values of n = (c'c), (x), and (p), to be 

, , 1 



Xes 



(54) 
(55) 



[l + ( K /2) 2 ] 

(P)oo = K«oo- (56) 



Thus the fixed point must satisfy: 



XES ( 1 



1 + (k/2) 2 \lRh L e 



4AS 00 ±4A<5 



(57) 



Following the derivation outlined in [Tol | , the stability of the fixed point can be found by looking at the eigenvalues 
derived from a linearised dynamics. From [Toj j. the system possesses one real eigenvalue and two complex conjugate 
eigenvalues which indicate the existence of a limit cycle. This limit cycle appears at a bifurcation of the fixed point 
which happens when the eigenvalues are purely imaginary. This condition appear at a critical value of xes'- 

(Al + A,K + (1) 2 + 1)A, K 
XES =Xh = tt (58) 

where A* = ^ L e~ 2A2 T2Ai + ^ Re 2\x±2\s^ 

Interestingly, the fixed point for the two spins are different and will result in a slightly different critical Xh values. 
However we choose to operate away from this regime in order to again avoid any possibility of noise upsetting the 
distinguishability of the resulting signals. 

A more straight forward method is to keep system such that it falls within the fixed point regime and proceed 
in calculating the difference in the steady state current. This is possible by a careful choice of values such that the 
coupling xes a re always below this critical values for both of the spins. This puts the system well in the supercritical 
regime. We can thus proceed by looking at the steady state at the fixed point and thus eliminating the possible error 
related to the existence of a stable limit cycle at the outer position of the doped fullcrcnc's oscillations. 

From (|54[) - l|5tj|) . we see that when A = 0, there is obviously no dependence of the current on the endohedral spin-state 
when Jr — Jl- This is not the case when the coupling is non- vanishing i.e. A > 0. This reflects the inherent asymmetry 
of the current transport from the source to the drain which is now mediated by an island which has, essentially, steady- 
state tunneling probabilities to either contact which are now effected by the endohedral spin dependent equilibrium 
position of the island. Taking, the steady-state semiclassical right tunnel current to be (ir)^ = 7r exp(±2<5A)(fi) 00 , 
we have 



■T \SC ^ +2A5 , ^ -2A<5 



i l R )§° l L e~ 2XS + 7i?e+ 2 ^ 



e 



(59) 



in the limit of large tr/tl (f r/Tl 3> e +4XS ), and large A<5. Thus we might expect an exponential dependance of the 
resulting current ratio on the endohedral spin-state. 

We can obtain the steady-state solution to Q51[l. by converting the resulting Master equation into a 4A^ 2 x 4A^ 2 , 
(where N is the size of the truncation of the Fock-state representation of the vibronic Hilbcrt space), complex 



9 



matrix and searching for the eigenstate (pm), with vanishing eigenvalue. We use the inverse power method to find 
the eigenvector with smallest eigenvalue In the following we mostly choose TV = 36, giving matrices of size 

4096 x 4096 ~ 256MBytes in memory. As the memory resources increase with TV 4 , we cannot realistically go beyond 
TV ~ 60. However, if the resulting large matrix is sparse then routines which only use matrix vector products can 
be used to determine the steady-state solution [l^. In Figure we see that although the matrix has many small 
elements, it is unclear that it is highly sparse. As the current (|53|1 . is exponential in x, the expectation value (in), 
could depend on very small matrix elements. To ensure convergence we will work with the full matrix in what follows 
and choose to truncate the bosonic vibrational Hilbert space at TV — 16, or TV = 32. 

In Figure^] we plot the "semiclassical" evolution of (J5UJ, and quantum steady-state expectation values when xes = 
3, K = 1, 7l = .1, 7# = .9, and A = S = .4. We see that the different endohedral spin states lead to very 
different phase space dynamics, differing island populations (n), and currents {iRjoa- One further notices a remarked 
difference between the magnitudes of the quantum and semiclassical expectation values for the currents in Fig. |j3Jc). 
This difference grows with the size of xes, and is primarily due to the extended nature of the quantum state in the 
vibronic phase space. With the availability of the full quantum steady-state solution we are also able to visualise the 
quantum state via quasi-distribution functions and in Figure we see that the quantum state is localised near the 
phase space origin. Moreover, from the compactness of the resulting Wigner or Husimi representations we can be 
relatively certain that the quantum steady-state solution is well represented with a Fock state truncation level set at 
TV = 32. As xes is increased towards the more physically realistic regime of xes ~ 20, one is justifiably concerned 
whether a truncation level of TV ~ 32, would remain sufficient to accurately represent the resulting quantum steady- 
state solution. We will see below however, that TV ~ 32, remains fairly accurate when xes = 10- In the section below 
we are primarily interested in knowing with relative certainty that our quantum simulations are highly accurate. In 
Figure IXll we superimpose the semiclassical phase space evolution on to the Wiger function of the quantum steady- 
state solution for the two cases of endohedral spin up and down. We note that due to a slight assymetry in the Wigner 
function the quantum expectation value for the steady-state equilibrium solution is not located at the peak of the 
Wigner function. 

We can see from this analysis that there is a significant difference between and (i R )- To characterize this 

difference more systematically we numerically evaluate the ratio (i R )/(i R ), as a function of tx and tr. We choose 
Xes = 3, R = 2, A = S = .4, and range through values of tx and jr, computing the resulting quantum steady-state 
current ratios for the two cases of endohedral spin up and spin down. The results are shown in Figure IXll and in the 
limit of large tr/tl (f r/T l 3> e +4XS ), the current ratio appears to asymptote, as it should to exp(4A£) ~ 1.9. By 
repeating these simulations at the reduced truncation TV = 16, and finding that the resulting current ratios are within 
10 -9 of those computed at TV = 32, indicates that our quantum steady-state simulations are extremely precise. 



C. Quantum Trajectories 

As mentioned above, going much beyond a Fock state truncation of TV ~ 60, with the full matrix representation of 
the Liouvillian is not possible without exotic computational resources. Another method of simulatin g th e dynamical 
evolution of the full quantum master equation (|51|) , is to use the method of quantum trajectories |lllll2l Il3l | . Briefly, 
in this method one represents the quantum master equation l|51|) , as an average of a stochastic Schroedinger equation 
over some stochastic noise. The temporal evolution of an initial state determined by this stochastic Schroedinger 
equation is known as a quantum trajectory. One can determine the time evolution of quantum expectation values 
of an observable by taking the expectation values of that observable within a single quantum trajectory and then 
averaging the expectation values found over numerous stochastic repetitions of the quantum trajectory. Since the 
method primarily involves the Schroedinger evolution of pure quantum states, the quantum evolution of an individual 
quantum trajectory can be simulated with much higher Fock state truncations TV, than in the above steady-state 
method. The down side is that to obtain precise expectation values one must average over a large number of repetitions 
of the quantum trajectory. 

The evolution of the stochastic Schrodingcr equation corresponding to l)51[). when TV = 0, is conditioned on the 
appearance of random quantum jumps. The information gained by the environment, in a sense, yields information 
which conditions the quantum state following the quantum jump. In l|51|) . there are three types of quantum jump: 
(I) — > c'El\iP) (unnormalised), which corresponds to an electron tunneling onto the island from the Fermi bath of 
the left electrode with some associated monitoring (via El), of the island's position state by the vibrational bath, (II) 
— > cEr\tJj), corresponding to an electron moving off the island into the Fermi bath of the right electrode with some 
associated monitoring (via Er), of the island's position state by the vibrational bath, and (III) — > a\ip), where a 
motional quanta escape's from the island's position state into the vibrational bath. In Figure HO we show an example 
of a single quantum trajectory of the coupled electronic- vibronic shuttle system corresponding to (|51|l . We choose no 
motional damping and set most of the parameters to vanish so we can illustrate the effects of the quantum jumps 
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clearly. At time points (A) in Figure a single electron jumps onto the island. Since this is more likely for negative 
5, due to the exponential dependance of El on position, we see a jump in the island's position towards the source 
contact located in the — x region. At time points (B), this electron now jumps off the island into the drain. Since 
again this is more likely closer to the drain due to the form of Er, we see a quantum jump in the island's position to 
larger values of +x. While the island is either empty/occupied, the non-Hermitian evolution suffered by the quantum 
trajectory causes the system to move towards the situation where the occupied/empty state becomes more and more 
probable. As a consequence of the exponential dependances of El and Er, this means that the oscillations grow 
towards larger +x, if the island is occupied, while they grow towards —x, if the island is empty. The net effect is to 
pump the vibrational motion of the shuttle during the single trajectory and even though this growth is interspersed 
with dips due to the movement on or off of an electron. Once averaged over many quantum trajectories, the resulting 
effect is to pump the average energy of the oscillation in an unbounded manner (if there is no motional damping) . 

We now consider the full-blown quantum trajectory simulation of the system using the parameter values 
Xes = 10, k = 2, 7l = .1, Jr = .8, and A = S = .4. We now have chosen a far larger value for the electrostatic 
force and are able to expand the Fock basis truncation to N = 100 using quantum trajectories. To determine the 
steady state quantum solution using a full matrix representation would now require over 26GBytes of memory! 
The resulting simulation averages over 2000 quantum trajectories and required many days to execute on a 2GHz 
Operton Dual processor 64-bit computer. In Figure we plot the right tunnel current (ir), as computed via 
quantum trajectories for the two endohedral spin states. We also plot the "semiclassical" dynamical evolution and 
the steady-state quantum expectation values for the currents where the latter is computed at the necessarily lower 
Fock truncation of N = 40. We can clearly see that the (incredibly computationally expensive), quantum trajectory 
simulation remarkably agrees very well with the values obtain from the steady-state numerics. It is not possible 
to plot out the Wigner function from the quantum trajectory simulation as this method can only yield the values 
of a small number of expectation values without running into the difficulty of requiring huge amounts of storage. 
However, the results shown in Figure ^\ strongly indicates that the technique of directly determining the quantum 
steady-state via the eigenvector of the Liouvillian matrix with vanishing eigenvalue, can extend to very large values 
of the electrostatic driving force xeSi without any significant lack of precision. As the physically relevant regime for 
Xes ~ 20, the behaviours we obtain here with xes ~ 3 — 10, should therefore also hold in the physically relevant 
regime. 



D. Estimating the Steady-State Quantum Current Noise 

There are a number of methods to deduce the noise power spectral density S(u>). We will use the methods outlined 
in [l|. 

We first represent the instantaneous right tunnel current as iR{t) = edN/dt, where dN(t), is a classical point process 
which represents the number (either zero or one), of tunneling events seen in an infinitesmal time, dt, and e is the 
electric charge. We can consider dN(t) to be the increment in the number of electrons N(t), in the drain electrode 
in a time dt. The steady-state current is computed by = E[i(i)]oo, where E[-], represents an average over the 
stochastic process describing the tunneling events. In our case this is given by, 

1% =lR(J[ce+^} Poc ) , (60) 

where the jump superoperator J'[A]B = ABA^. The fluctuations in the current are quantified by the two-time 
correlation function 

G(t) = E[i R (t)i R (t + t) - i^U = ei%6(t) + (i R (t)i R (t + t))%° . (61) 
From the theory of open quantum systems [T3| . one can show that 

G(t) = ei%S(t) + e 2 I Tr[J R e ct jR Poo } - Tr[J RPoo } 2 \_ , (62) 

where Jr = t/jr c cxp(— Ax) (see l|21[l). and the Liouvilian evolution, exp(£i), is according to the master equation 
(|4*T)l . In the case of a decoupled system (A = 0), this leads to a normally ordered correlation function while one 
obtains an antinormally ordered correlation function if one were studying the noise of the left tunnel current [T^ . For 
the coupled system the interpretation of the operator ordering is more complicated. 
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The current noise spectral density is given by 



S(uj) = 2 dtG(t) e iut . (63) 

J — CO 

In terms of our rescaled variables and parameters, r = cuot, and (|46|1 . we set f — toot, to obtain, 

( _ _ A OC 

G(f) = ei%8(f )lj + e 2 cu 2 j 2 R { Tr[J[ce+ x *]e ct J[+ce +X£ } Poo ] - Tr[J[ce +Xs ] Poo ] 2 \ , (64) 
where now the Liouvillian evolution is given by l|51[) . The spectral density is now given by 

e i{^ )f d ~ f T T [j[ ce +^] e ^ j[ ce +^] Poo ] - Tc[J[ce +Xs ] Poo ] 2 } . (65) 
From this we can derive the Fano factor, F = S(ui)/ (2ei R ), to be 

F(o;/a;o) = 1 + 4^ e *(-M>)f df {Tr[ l 7[ce+^]e £ *J[ C e +As ] Poo ] - Tr{J{ce+ Xs } Poo } 2 } . (66) 

The Fano factor F gives information on the statistics of the tunnel current noise. If F = 1 then the noise is 
completely uncorrelated/Poissonian, i.e. white noise. If F > 1, then the noise is known as super-Poissonian and 
the tunnel events are bunched. If F < 1, then the noise is known as sub-Poissonian and the tunneling events are 
anti-bunched, i.e. there is little chance that a second tunnel event will closely follow a previous event. Quantum 
correlations within the coupled system are primarily responsible for F ^ 0. In the case of the decoupled system, 
In Figure [XTI we compare the results for the decoupled system (A = 0), with the standard results for a two-state 
sequential tunneling process , 

7 oo _ 1L1R p _ a + ll ^ 



1L+1R 1 (7£ + in) 2 

In Figure [TH we display the current noise spectra for various frequencies, i.e. 

*VM0 - 1 - ■ (68) 

and see that typically the noise is sub-Poissonian. In Figure IT21 wc contour plot the DC (or zero frequency), Fano 
factors for the decoupled and coupled system and observe that the coupled system displays stronger anti-bunching 
when ji ~ 7a, than in the decoupled case. However in the parameter region where we are interested, jr 3> 7l, the 
DC Fano factors of the decoupled and coupled system are very similar. 



E. Measurement Time 



cndohedral spin-states, is the measurement time r m , required to differentiate the two currents i R °° and i R °°- This is 



The final quantity we must determine, given the DC current signals and current noise associated with the two 
idohedral spin-states, is 1 
given by the expression ^ 



2(Aiy 



(69) 



where Si, and S2, are the DC current noise spectral densities for the two cases of spin up/down while A J = \ (i R { 



T 00 \ 



1 00 \ 

R 



. In Figure H31 wc plot log 10 r m , for various jl, Jr- We see clearly that when jr ~ 7^, r m must be very large 
to discriminate between the very small difference in steady-state currents. In the parameter regions we are concerned 
with, i.e. 7r/tl S> 1, this measurement time is r m ~ 10 1 — 10 2 , in the natural units of the oscillator (cj = 27r x 1GHz), 
or a few tens of nanoseconds. Thus in this parameter region, the measurement time is very realistic. 
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X. SUMMARY 

In summary, we have shown that the right-hand tunnel current of a strongly damped single electron shuttle operating 
in the Coulomb blockade regime has very high sensitivity to small equilibrium position displacements of the shuttle. 
We considered a N@Cgo molecular shuttle whose endohedral spin can exert a force on the molecule when placed in 
a large magnetic field gradient. By placing a large nanoscopic magnet nearby, a tiny ~ ±10~ 15 N force is suffered 
by the molecular island arising from the endohedral spin state S z = ±3/2. This small force alters the equilibrium 
position of the island in the presence of the electrostatic driving force, motional damping, and harmonic restoring 
force. Through the extremely sensitive (exponential), dependence of the conductance on the tunnel separations, the 
shift in the equilibrium position of the island profoundly influences the current flow. We determined the current noise 
spectral density and the measurement time required to distinguish between the two DC steady-state current signals 
(spin up and down), in the presence of steady-state quantum noise. 

The results we have found would strongly indicate that there are realistic parameter regimes where the spin state 
dependent currents are distinguishable within several tens of nanoseconds. Thus this device should thus be capable 
of spin-detection. Alternatively, similar devices (without an endohedral spin), should be capable of discriminating 
small, ~ ±10 _15 N , static forces on a Cqq island. 
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FIG. 1: Schematic of the N@C6omolecular shuttle. NM labels the nanomagnet, while the left and right tunnel contacts are 
held a the voltages Vl and Vr respectively. 




FIG. 2: Schematic of the nanomagnet, a small spherical ferromagnetic particle NM, of radius R, a distance d from the shuttle 
S. 




FIG. 3: Sparsity of the Liouvillian superoperator. Plot of log 10 of \pij\, when expanded into a AN 2 x 2N 2 , N = 16 matrix. 
The steady-state solution is the eigenvector with vanishing eigenvalue of this matrix. 
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FIG. 4: Semiclassical and steady-state quantum expectation values for the case \es = 3, R = 1, 7l = .1, 7h = -9, and 
X = 5 = A. (A) semiclassical phase space evolution (in displaced picture), (a)/(b) correspond to spin up/down or ±5; (B) 
graphs of (c'c), for the two endohedral spin states (a) and (b), and where solid curves are the semiclassical solutions and the 
dashed are the quantum steady-state expectation values with N — 32; (C) graphs of (is), for the two endohedral spin states 
(a) and (b), and where solid curves are the semiclassical solutions and the dashed are the quantum steady-state expectation 
values with N = 32. 
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FIG. 6: Semiclassical dynamics overlayed with the quantum steady-state. Graphs (A) and (C) are the semiclassical (blue/black 
curves), quantum steady-state phase point ((i)oo, (p)oo), as a red symbol, and the contours of the quantum steady-state Wigner 
function, for the case of endohedral spin up and down. The blue curve is the appropriate semiclassical phase-space trajectory 
for the spin up(down), while the black curve is the other case shown for comparison. Graphs (B) and (D) are same except the 
Husimi function is used. 



17 




FIG. 7: Graph of (i R )cc/{i R )co, computed via the quantum steady-state for a range of values of the left and right tunnel 
amplitudes, jl and tr. For large 7_r/7l, the ratio tends to exp(+4A5). 




FIG. 8: Quantum trajectories simulation (1 trajectory with a basis of 100 Fock states), of (I41L for the case when = jl ~ 1, 
A = S — 0.4, R, = n — 8 = rj = \es = 0. Graph (a) (x/x ) = (x) vs. r, (b) (c'c), (c) (a^a). See text for description of quantum 
jumps at times A and B. 
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FIG. 9: Graphs of the right tunnel current given as a function of rescaled time t, for the two spin-states (±5). Graphs 
(a-b)/ (blue-red) is (in) computed via a quantum trajectories solution of J4H . for spin up/down (2000 trajectories). Graphs 
(c-d) (turquoise-orange) are the re-plotted "semiclassical" solutions for the two cases of the shuttle spin state. Graphs (e- 
f) (magenta-purple) are {in} computed via the steady-state solution for 1411 . on a N — 40, Fock basis representation. Again we 
have set R = 2, \es = 10, 7l = .1, "fR = -8, A = 5 = .4. 
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FIG. 10: Graph of Fano function F, for the decoupled system A = 0, k 
^ = (7i+7l)/(7i+7fl) 2 - 



= .1, Xes = 0. Numerical results agree with 
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FIG. 12: Contour plots of the DC Fano factor for the (a) decoupled system with \es = 0, k = .1, and (b) the coupled system 
with A = .4, xes = 3, R — 3. 
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FIG. 13: Contour plot of the log 10 (r m ), of the measurement time from Equation 1691 . required to distinguish the two spin 
states from the current for the coupled system with A = 5 = .4, \es = 3, R = 3. 



